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Abstract 

In this paper, we develop a multi-group epidemic framework via virtual dispersal where the 
risk of infection is a function of the residence time and local environmental risk. This novel 
approach eliminates the need to define and measure contact rates that are used in the tradi¬ 
tional multi-group epidemic models with heterogeneous mixing. We apply this approach to a 
general n-patch SIS model whose basic reproduction number IZo is computed as a function of a 
patch residence-times matrix P. Our analysis implies that the resulting n-patch SIS model has 
robust dynamics when patches are strongly connected: there is a globally stable endemic equi¬ 
librium when TZo > 1 while the disease free equilibrium is globally stable when TZo < 1. Our 
further analysis indicates that the dispersal behavior described by the residence-times matrix 
P has profound effects on the disease dynamics at the single patch level with consequences 
that proper dispersal behavior along with the local environmental risk can either promote or 
eliminate the endemic in particular patches. Our work highlights the impact of residence times 
matrix if the patches are not strongly connected. Our framework can be generalized in other 
endemic and disease outbreak models. As an illustration, we apply our framework to a two- 
patch SIR single outbreak epidemic model where the process of disease invasion is connected 
to the final epidemic size relationship. We also explore the impact of disease prevalence driven 
decision using a phenomenological modeling approach in order to contrast the role of constant 
versus state dependent P on disease dynamics. 


1. Introduction 

Sir Ronald Ross must be considered the founder of mathematical epidemiology [52] despite the fact 
that Daniel Bernouilli ( 1700 - 1782 ), was most likely the first researcher to introduce the use of math¬ 
ematical models in the study of epidemic outbreaks [S, 28] nearly 150 years earlier. Ross’ appendix to 
his 1911 paper [ ] not only introduces a nonlinear system of differential equations aimed at captur¬ 

ing the overall dynamics of malaria contagion, a disease driven by the interactions of hosts, vectors 
and the life-history of Plasmodium falciparum, but also includes a tribute to mathematics through 
his observation that this framework, his model, may also be used to model the dynamics of sexually 
transmitted diseases [52]. Ross’ observation has motivated the use of mathematics in the study of 
the impact of human social interaction on disease dynamics [9, 1 , 21 , 23, 22, 3 , 16, 41, 42, 43, 58]. 

In fact, Ross’ work introduced the type of frameworks needed to capture and modify the dynam¬ 
ics of epidemic outbreaks; new landscapes where public policies could be tried and tested without 
harming anybody, complementing and expanding the role that statistics plays in epidemiology. 
Suddenly scientists and public health experts had a “laboratory” for assessing the impact of trans¬ 
mission mechanisms; evaluating, a priori, efforts aimed at mitigating or eliminating the deleterious 
impact of disease dynamics. 


AMS 2000 subject classifications: Primary 34D23, 92D25, 60K35 

Keywords and phrases: Epidemiology; SIS-SIR Models; Dispersal; Residence Times; Global Stability; Adaptive 
Behavior; Final Size Relationship. 


1 



2 Derdei Bichara, Yun Kang, Carlos Castillo-Chavez, Richard Horan and Charles Perrings 

The study of the dynamics of communicable disease in metapopulation, multi-group or age- 
structure models has also benefitted from the work of Ross. Contact matrices have been used in 
the study of disease dynamics to accommodate or capture the dynamics of heterogeneous mixing 
populations [1, 19, 29, 35]. The spread of communicable diseases like measles, chicken pox or rubella 
is intimately connected to the the concept of contact, “effective” contact or “effective” per capita 
contact rate [25, 35]; a clear measurable concept in, for example, the context of sexually transmitted 
diseases (STDs) or vector-borne diseases. The values used to define a contact matrix emerge from 
the a priori belief that contacts can be clearly defined and measured in any context. Their use in 
the context of communicable diseases is based often on relative rankings; the result of observational 
subjective measures of contact or activity levels. For example, since children are believed to have 
the most contacts per unit of time, their observed activity levels are routinely used to set a relative 
contact or activity scale. Traditionally, since school children are assumed to be the most active, 
they are used to set the scale with the rest of the age-specific contact matrix usually completed 
under the assumption of proportionate (weighted random) mixing (albeit other forms of mixing 
are possible [2, 9, 1 , 20, 29, 3 ] and references therein). In short, mixing or contact matrices are 
used to collect re-scaled estimated levels of activity among interacting subgroups or age-classes; a 
phenomenological estimation process based on observational studies, surveys, and various appealing 
definitions of contact [51 ]. Our belief that contact rates cannot, in general, be measured in satisfac¬ 
tory ways for diseases like influenza, measles or tuberculosis, arises from the difficulty of assessing 
the average number of contacts per unit of time of children in a school bus, or the average number 
of contacts per unit of time that children and adults have with each other in a classroom or at the 
library, per unit of time. In some cities in Latin America, some individuals spend 2-4 hours per 
day as users of mass transportation systems, some traveling in packed subway cars or as “sardines” 
in small buses. How packed these modes of transportation are as a function of time of the day or 
day of the week can be observed but has not been uniformly quantified in terms of contacts (or 
age-specific contacts) per unit of time by different observers. The issue is further confounded by 
our inability to assess what an effective contact is: a definition that may have to be tied in to the 
density of floating virus particles, air circulation patterns, or whether or not contaminated surfaces 
are touched by susceptible individuals. In short, defining and measuring a contact or an effective 
contact, turns out to be incredibly challenging [5( ]. That said, experimental methods may be used 
to estimate the average risk of acquiring, for example, tuberculosis (TB) or influenza, to individuals 
that spend on the average 3 hours per day in public transportation, in Mexico City or New York City. 


In this paper we propose the use of residence times in heterogeneous environments, as a proxy 
for “effective” contacts over an u x” windows in time. Catching a communicable disease would of 
course depend on the presence of infected/infectious individuals (a necessary condition), the level 
of “risk” within a given “patch” (crowded bars, airports, schools, work places, etc), and the time 
spent in such environment. Risk of infection is assumed to be a function of the time spent in pre¬ 
specified environments; risk that may be experimentally measured. We argue that characterizing 
a landscape as a collection of patches defined by risk (public transportation, schools, malls, work 
place, homes, etc) is possible, especially if the risk of infection in such “local” environments is in 
addition a function of residence times and disease levels. Ranking patch-dependent risks of infection 
via the values of the transmission rate (/3) per unit of time, may therefore be possible and useful. 
The reinterpretation of (3 and the use of residence times move us away from the world of models 
that account for transmission via the use of differential susceptibility to the world where infection 
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depends on local environmental risk. 

Consequently, we introduce a residence times framework in the context of a multi-group system 
defined by patch-dependent risk (defined by /3). We study the role of patch residence times on dis¬ 
ease dynamics within endemic and single outbreak multi-group scenarios. Specifically, the study of 
the impact of patch residence times (modeled by a matrix of constants) on disease dynamics within 
a Susceptible-Infected-Susceptible (SIS) framework is carried out first, under the philosophy found 
in [ , 11, 13, 15, 17, 21, 3 , 14]. Individuals move across patches as a function of their assessment 

of relative levels of infection in each area (studies using alternative classical approaches are found in 
[15, 16, 32, 56]). Generalizations are explored through simulations of the two-patch SIS model with 
state-dependent residence times within our framework. The results are compared to the disease 
dynamics generated by constant residence times. 

The rest of this paper is organized as follows: Section 2 introduces a general n patch SIS model 
that accounts for residence times. Theoretical results on the role of residence times matrix (P) on 
disease dynamics are carried out using the residence times dependent basic reproduction number 
Ho (P). Patch-specific reproduction numbers 7£g(P), i = 1,..., n are used to highlight the impact of 
residence times matrices on cases that includes non-strongly connected patch configuration. Section 
3 explores, through simulations, the dynamics of the SIS model under a state-dependent residence 
times matrix in a two-patch system; P = P(/i,/ 2 )- That is, when the decisions to spend time in 
a patch are a function of patch-disease prevalence. Section 4 highlights our framework in the case 
of a two patch single outbreak SIR model following the work of Brauer [10, 17], and discusses the 
role of P on the final epidemic size. Section 5 collects our observations, conclusions and discusses 
future work. The detailed proofs of our theoretical results are provided in the Appendix. 


2. A general n-patch SIS model with residence times 

A general n-patch SIS model with residence time matrix P is derived. The global analysis of the 
model is carried out via the basic reproduction number TZq. We also include patch-dependent dis¬ 
ease persistence conditions. 

2.1. Model derivation We model disease dynamics within an environment defined by n patches 
(or risk areas) and so, we let Ni(t),i = 1,2 ...,n denote resident population at Patch i at time t. 
We assume that Patch i residents spend pij E [0,1] time in Patch j, with Ylj=iPij = f° r each 
i = 1,..., n. In extreme cases, for examples, we may have, for p^ = 0, i 7 ^ j, that is Patch i residents 
spend no time in Patch j while Ylj^iPij = 1 (° r equivalently p l% = 0) would imply that Patch i 
residents spend all their time in Patch j (with j = 1 ,,n and j 7 ^ i ) even though their patch is 
(labelled) i. In the absence of disease dynamics, the population of Patch i residents is modeled by 
the following equation: 

= bi- diNi (1) 

where the parameters b {, di represent the birth rate, and the natural per capita death rate in Patch 
i , respectively. Hence, the Patch i resident population approaches the constant ^ as t—> 00 . 
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In the presence of disease, we assume that disease dynamics are captured by an SIS model, 
thus, the Patch i resident population is divided into susceptible and infected classes, represented by 
Si, It , respectively, with Si + /, = N,. We further assume that (a) there is no additional death due 
to disease; (b) the Patch i Infected resident population recovers and goes back to the susceptible 
class at the per capita rate 7 f, (c) the residence time matrix P = collects the proportion 

of times spent by z-residents in j -environments, i = 1,..., n and j = 1,..., n. The disease dynamics 
are therefore described by the following equations: 

Si = bi — diSi + 7 ili - J2 T j=i(Si infected in Patch j) 

Ii = Y^=i( s i infected in Patch j) - 7 *U - d,I t (2) 

Ni = bi-diNi. 

We model Si infection within Patch j in the following way: 

• Since each pjj entry of P denotes the proportion of time that Patch i residents spent mingling 
in Patch j, we have that: 

— There are PijNi = pijSi + Pijh Patch i residents in Patch j on the average at time t. 

— The total Patch j , the total effective population is Pkj^k- of which Ylk=iPkjIk are 

V-m r 

infected. Hence, the proportion of infected individuals in Patch j is j2™ =1 p * j .N k an< 4 we ^ 
defined, as long as there exists a k such that pkj > 0 ; so that the population in Patch j 
is nonzero. 

• Hence, the [Si infected per unit of time in Patch j] can be represented as the product of the 
following three items: 



the risk of infection in Patch j 


X 



Susceptible from Patch i who are currently in Patch j 


1 Pkj Ik 

ELi PkjNk 

'-V-' 

Proportion of infected in Patch j 

The transmission takes on a modified frequency-dependent form that depends on how much 
time individuals of each epidemiological class spend in a particular area, and where fi 3 differs 
by patch to reflect spatial differences in potential infectivity. More precisely, ,6j is assumed 
to be a patch-specific measure of disease risk per unit of time with its effectiveness tied in to 
local environmental and sanitary conditions. Therefore, 

[Si infected per unit of time in Patch j] = f3j x p^Si x J ~^ =1 - — ^ (3) 

l^k= 1 PkjNk 


provided that there exists k such that pkj > 0. 
Model (2) can be rewritten as follows: 


Si = bi - diSi + 7 ili - E" = 1 (ft x PijSi X 
ii = Sj=l {Pj X PijSi x ) _ Ti-fi — dili, 

Ni = bi - diNi, 


( 4 ) 
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with, the dynamics of the Patch i resident total population modeled by the equation: Ni(t) = 
bi~ djiVj(t), where 5) + = A T t , which implies that IVj(i)—as t—> + oo. Theory of asymptotically 

Cti 

autonomous systems for triangular systems [24, 57] guaranties that System (4) is asymptotically 
equivalent to: 


h — I PjPij ( d, I' 


T,k=l Pkjlk 

T,k=l Pkj J7 


— ( 7 * + di)Ii 


■>n _ PjPi 

~vn 

-*k 


iPjj 1 I (h. _ t\ V™ P jPij ^k=l,k^iPkjIk _ / , \ T 

h )^j =1 w- v ..hL Vh + hJJ-i 

lPkjd kJ ' ' I^k=lPkj dk 


(5) 


= IS* 

' ' y ^t=i«3 d fc j 

for z = 1 , 2 ,..., n, with residence times matrix P = (pij)jZ j’"’” satisfying the conditions: 
HP1. At least one entry in each column of P is strictly positive; and 
HP2. The sum of all entries in each row is one, i.e., YYYj =1 Pv = 1 f° r all *■ 


2.2. Equilibria, basic reproduction number and global analysis To analyze the system, 
we investigate the basic reproduction number of the system with fixed residence times to better 
understand its properties in the absence of behavioral responses to risk. We let B = (j3i,/3 2 , • • •, /3 n ) 
define the risk of infection vector; (3i a measure of the risk per susceptible per unit of time while in 
residence in Patch i. 


Letting S = (Si, S 2 ,..., S n Y , I = (h,h, • • •, I n f , N = 

( ELi \ 


in h 

d 1 , d 2 


1 7 ? * * * 5 


dr, 


, and TV = P^TV = 


E n bi, 

k=l 


V Pknd k ) 


. Then System (5) can be rewritten in the following compact (vectorial) form: 
i = diag(lV — I)Pdiag(S)diag(iV)^ 1 P f / — diag(d/ + 7 j)I ( 6 ) 


with state space in R”. System ( 6 ) has the compact set Q = {I > Ok™, I < N} as its global 
attractor. This implies that the populations involved are “biologically” well-defined since solutions 
of ( 6 ) will converge to and stay in fi. We therefore restrict the dynamics of ( 6 ) to the compact set Q. 


The analysis of System ( 6 ) is naturally tied in to the basic reproductive number TZq [ !7, 55]; the 
average number of secondary cases produced by an infected individual during its infectious period 
while interacting with a purely susceptible population. IZo is given by (see the detailed formulation 
in Appendix): 

IZq = p(— diag(lV)Pdiag(S)diag(iV) _ 1 P i l/ _1 ) 
where V = -diag (d/ + 7 /), dj = (d 1} d 2 , ■ ■ ■, d n Y and 7 / = ( 71 , 72 , • • •, 7 nY■ 


( 7 ) 
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The basic reproduction number IZq is used to establish global properties of System (6). For the 
relevant literature on global stability for multi-group or metapopulation models, see [5, 46, 47, 48, 54] 
and the references therein. We define the disease free equilibrium (DFE) of System (6) as I* = 0®™ 
and the endemic equilibrium (when 1Z o > 1) as I where all components are positive. By using the 
same approach as in [ 16, 48], we arrive at the following theorem regarding the global dynamics of 
Model (6). 

Theorem 2.1. [Global dynamics of Model (6)] Suppose that the residence times matrix P is 
irreducible, then the following statements hold: 

• IflZ o < 1, the DFE I* = O^n is globally asymptotically stable. IflZo > 1 the DFE is unstable. 

• IflZ o > 1, there exists a unique endemic equilibrium I which is GAS. 

Remarks: The detailed proof of Theorem 2.1 is provided in Appendix B. These results imply that 
System (6) is robust, that is, disease outcomes are completely determined by whether or not the 
reproduction number IZq is greater or less than one. The results of Theorem 2.1 while powerful, 
do not provide easily accessible insights on the impact of the residence matrix P on the levels of 
infection within each patch. 


Direct insights on the effects of P, are derived by focusing on the levels of endemicity within each 
patch. The following two definitions help set the stage for the discussion: 

• The basic reproduction number for Patch i in the absence of movement (pa = 1 or Pij = 
0), SIS model, is defined as 1Z q = d .+ ■, which determines whether or not the disease will be 
endemic in Patch i. In short disease will die out if TZq < 1 with a unique endemic equilibrium, 
that is GAS, if IZq > 1. 

• The basic reproduction number associated with Patch i, under the presence of multi-patch 
residents, is defined as follows: 






. b k 
Ft 


577=1 hjPii 


v ■ bi 




di+'Yi 


~ K X YTj= 1 ( ft ) Pij 




di+'n 


£ 7 = 1^2 


We explore the role that IZq (P) plays in determining the impact of all residents on disease dynamics 
persistence in Patch i in the following theorem. 

Theorem 2.2. [The endemicity of disease in Patch i] Assume that the residence times matrix 
P satisfies Condition HP1 and HP2 but that some of its entries can be zeros. 

• IflZ q(P) > 1, then the disease persists in Patch i. 

• If the following conditions hold: 


H: pkj = 0 for all k = 1,.., n, and k i, whenever pij > 0, 


ftj(P) = Kx'E 

3 =1 




n 


= ^o*E 

3 = 1 



Pij- 


then we have 
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Thus, when Condition H holds and TZ l 0 X Ybj =1 Pij < 1, then endemic levels of disease 
cannot be supported in Patch i. That is, 


lim Ii(t) = 0. 
t—>oo 


Remarks: The detailed proof of Theorem 2.2 is provided in Appendix C. The results of Theorem 
2.2 give insights on the role that the infection risk (measured by B) and the residence time matrix 
(P) have in promoting or suppressing infection. Further, a closer look at the expression of the 
general basic reproduction number in Patch i, namely 


ni i 


= K 


X > I — 


E 

3 = 1 


Pij 


Pvd, 




leads to the following observations: 

1. The movement between patches, modeled via residence time matrix P, can promote endemic- 

ity: For example, if 1Z 1 0 = < 1, i.e., there is no endemic disease in patch i. Then, the 

presence of movement connecting Patch i to possibly all other patches can support endemic 
disease levels in the following ways: 

• Via the presence of high risk patches, that is, there exists a patch j such that is large 
enough. For example, letting pki = 1 /n for all k , l with the total population in each 

v - "' n n 

patch being the same = K for all k ; K a constant) then 77 .q(P) = J and 

consequently, if J2j=i Pj > then Patch i will promote the disease at endemic levels. 

• Whenever individuals spend more time in high risk than in low risk patches. For example, 

in the extreme case, p^ = 1 with we have that 7 £q(P) > 1 , and thus, endemic 

disease levels in Patch i can be supported. Patch j (j = 1,... ,n and j 7 ^ i) can therefore 
be considered the source and Patch i (i / j) the sink [3, 4, 5, 6 , 47, 48, 54, 53]. 

2 . Under the assumption TZ' l 0 > 1, for an isolated Patch i, conditions that lead to disease extinc¬ 
tion in the same Patch i under the movement can be identified. According to Theorem 2 . 2 , 
Condition J should be satisfied and so the expression of 77 .q(P) reduces to 


t ^( p ) = n l 0 xJ2 

3 = 1 



( J, >! 1 j 


n 


— K-b x E 

3 = 1 



Pij. 


Therefore, the only way to have the value of 77 .q(P) be less than one, would be when the 
amount of time spent in Patch i is such that Xq=i (/3 1 ) ^ i/J < E(p) • Therefore, we conclude 
that the synergy between the residence time matrix P and the existence of sufficient low risk 
patches (i.e., (3j <C fif) can suppress a disease outbreak in Patch i. 


3. Two patch models: state-dependent residence times matrix 

We now extend the analysis of disease dynamics to the case where susceptible individuals respond 
to variations in risk in an automatic way. In particular, we consider the case when susceptible 
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individuals make programmed responses to variations in disease risk, and do not choose their re¬ 
sponse to optimize an index of wellbeing ( see for example [10, 13, 15, 17]). While this may not 
be a very good approximation of disease risk management in real systems, it enables us to explore 
the implications of certain types of phenomenologically modeled behavioral responses by assuming, 
for example, that the proportion of time spent in a particular patch depends on the numbers of 
infected individuals on that particular patch; that is P = P(/i,/ 2 ). 


Possible properties of the proportion of time spent by resident of Patch i into Patch j , i ^ j, ( pij ) 
may include: increases with respect to the growth of infected resident in patch i (/*), or decreases 
with respect to infected resident in patch j ( Ij ). Mathematically, we would have that 


Ij) 

Bln 


< 0 and 


dpij ( 1% 11j 

Oh. 


> 0 . 


In a two-patch system, the use of the relationship pij(I±,l 2 ) + Pji(h, h) = 1, reduces the above 
four conditions on P, to the following conditions: 


dpn(h,l2) , dp22(h,l2) 

< o and 1 \ T ; < 0. 
oI\ 0 I 2 

Examples of functions Pij{I\ , I 2 ) with these properties include, 


and 


Pi 2 {h,h) = q- 12 .1 1 + 7l T and p 2 i(h, h) = 021 1 + ^ 

1 + il + 12 l 12 


n T\ °11 + Vllh + h J /T T\ °22 + I\ + V 22 I 2 

Pii{h,l 2 ) = — , , J , r - and P 22 {h,l 2 ) = 


1 + h + I 2 


1 + h + I 2 


2 

where &ij are such that Oij = 1 . 

3 = 1 

More complex behavioral response formulations may also depend on the states of total popula¬ 
tions Ni and N 2 , but the current specification captures important components of risk (infections) 
and allows us to retain the asymptotic equivalence property applied in the case of fixed residence 
times. Hence, using the same notation as in System ( 6 ) leads to the following two dimensional 
system with P = P(Ji, .Z 2 ): 


where 


il = X(Ii,l2)(^-Ii)h + Y(Ii,l2)(^-Ii)l2 - (di +71 )Iu 
h = Y(h,I 2 )( | - h)h + Z(h,I 2 ){% - h)h - (d 2 + 72K2, 


X(h,I 2 ) = 


Pipli(Ii,h) 


+ 


fhp^A- 12 ) 


Pnihih)^ +P2i{Ii,h)^ Pi2{h,l2)j; + P22(h, h)^; 


d2 


di 


d2 


Y(Il,I 2 ) 


/3lPll(h, l2)P2l{Il, I2) 

Pll(/l,/ 2 )|+P 2 l(/l,/ 2 )| 


P2Pl2(h,l2)p22(h,l2) 

Pl2(hJ2)^ +P22{Ilj2)^' 


( 8 ) 
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and 


Z(h,h) = 


Piph(h,h) 


+ 




Pll{hj2) b ±+P2l{h,h) b ± Pl2(h,h) b ± 1 +P22(.h,h) b £ 


I dl I yZL\*L,*ZJ d2 

where A'(/i,/ 2 ), T(7l,/ 2 ) and Z{I\,l 2 ) are positive functions of I\ and I 2 . 


The basic reproduction number 1Z 0 is the same as in the previous section since it is computed at 
the infection-free state, i.e. 

TZq = / o(diag(iV)Pdiag(S)diag(7V)" 1 P f (—f/ -1 )) 


where, in this case, we have that P = 


0"ll 

<712 

<721 

<722 


and aij=pij(0 ,0), V{i, j} = {1, 2}. 


The properties of positiveness and boundedness of trajectories of System (6) are preserved in 
System (8). In addition, System (8) has a unique DFE equilibrium whose local stability is deter¬ 
mined by the value of the (uncontrolled) TZq\ the DFE is locally asymptotically stable if TZq < 1 
while it is unstable if 1Z 0 > 1. 


Let us consider whether System (8) can have a boundary equilibrium such as (0,7 2 ) or (7i,0). 
The assumption that System (8) has such a boundary equilibrium (0,7 2 ) with 7 2 > 0 implies that 
r(0,/ 2 ) = 0. Since pn(0, 12 ) = and p 2 2(0,1 2 ) = 022 , we deduce that 


Y( 0 ,h) = 


1+/2 

/3icr 2 i(crii + I 2 ) 


+ 


/?2<712<722 


°-^ b -± + <?2l b ±(l + h) ' <U2^ + ^22^(1 + I2) 


1+h <h 1 

This indicates that Y (0, / 2 ) = 0 if and only if ct 2 i = 0 and o \2 = 0, which requires that: 


P 12 = P 21 = 0, and pii = P 22 = 1. 

A similar arguments can be applied to the boundary equilibrium (li,0). Therefore, we conclude 
that System (8) will have a boundary equilibrium f (0, J 2 ) or (7i,0)) only in the trivial case of 
isolated patches, that is, where there is no movement between two patches. This conclusion differs 
from the state-independent residence matrix model (6), since for example, the two-patch model (6), 
according to Theorem 2.1, boundary equilibrium (0,7 2 ) or (7i,0) can exist when pn = p 22 = 0 
(pi 2 = P 21 = !)• 


To illustrate the difference between the state-dependent residence matrix model (8) and the state- 
independent residence matrix model (6), we look at the situation when an = d 22 = 0, oq 2 = cr 2 1 = 1 
( P 11 = P 22 = 0,pi 2 = p 2 i = 0 for the state-independent residence matrix model (6)). Under the 
condition of an = cr 22 = 0,ui 2 = cr 2 i = 1, we have Model (8), that 


and 


Pl2(/l,/ 2 ) 


l+/l 
1 + 7i + J 2 


and P2i(7i,7 2 ) 


1 + 7 2 
1 + 7i + I 2 


h 


h 

1 + 7i + I 2 


Pn(h,h) 


1 + Ii + I 2 


and p 22 (/i,/ 2 ) 
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This difference has significant impact on disease dynamics (see Fig 1(a) and Fig 1(b), red curves). 

In Fig 1(b), we see that the infection in Patch 2 (high risk) persists in the state-dependent case 
whereas it dies out when P is constant. That is due to the fact that pa(Ii,l 2 ) will not equal zero 
whereas Pij ( I \ , I2 ) with i j may. For the constant residence times matrix, the dynamics of the 
disease in each patch is also independent, where people in patch i infect only susceptible in patch 
j with i ^ j. In Fig 1(b) (red solid curve), we observe that the disease dies out in Patch 2 with 
TZq = = 0.8571. For the state-dependent case, unless there is no disease in both patches or 

one disease-free Patch, the proportion of time residents spend in their own patch is nonzero. This 
leads the disease to persist in both patches if TZq > 1 (see Fig 1(b), red dashed curves). However, 
even in this case, the disease dies out in both patches if TZq < 1 (See Fig 3, red curves, for instance). 



(a) Dynamics of the disease in Patch 1. If there is no (b) Dynamics of the disease in Patch 2. In the high 
movement between the patches (blue curves), the dis- mobility case, the disease dies out (solid red curve) 
ease dies out in the low risk Patch 1 in both approaches for P constant, with TZg = 0.8571, and persists for P 
with TZq = 0.7636. state-dependent (dashed red curve). 


Figure 1 . Coupled Dynamics of I\ and I2 for constant pij (solid) and state dependent pij (dashed). The red lines is 
case of high mobility, i.e. P12 = P21 = cr 12 = 021 = 1 - The black lines represent the symmetric case, i.e: p 12 = P21 = 
<Ji2 = 021 = 0.5 and the blue line represent the polar case, i.e.: P12 = P21 = 012 = 021 = 0. 



Figure 2 . Coupled Dynamics of I1+I2 for constant pij (solid) and state dependent pij (dashed). The overall prevalence 
is higher if if the residence times is symmetric (solid and dashed black curves). The black curves represent the 
symmetric case ( P 12 = P 21 = <712 = 021 = 0.5 ), and the blue line represent the polar case ( P 12 = P 21 = <r 12 = 021 = 0) 
and red curves represent high mobility case (p 12 = P 21 = cr 12 = (721 = 1 ). 
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3.1. Applications and comparisons: the two patch cases The analytical results of the 
global dynamics on the asymptotic behavior of Model ( 8 ) are still unresolved. Hence, we ran simu¬ 
lations to gain some insights on the role of P(/i, I2) on endemic dynamics. We observe that trajec¬ 
tories converge towards an endemic equilibrium whenever IZo > 1 ; however, there are substantial 
differences in the transient dynamics generated by state-dependent when compared to 

those generated with a constant residence times matrix. 

Unless stated otherwise, we suppose the following generic values for the simulations: /?i = 
0.3, P 2 = 1-2, b\ = 9 , di = 1/7, &2 = 9 , d ,2 = 1/10 and 71 = 72 = 1/4. From a selected of 
simulations, it is observed that: 

1. For the symmetric case where p \2 = P 21 = 0.5, the disease is endemic in both patches as pre¬ 
dicted by Theorem 2.1 since IZo = 2.0466. For the state-dependent case, simulations suggest 
(Fig 1(a) and Fig 1(b), black dashed curves) that trajectories tend to be endemic in both 
patches. However, the level of endemicity is lower than the constant case in Patch 1 (low risk 
patch) and is greater in Patch 2 (high risk patch). 

2. Fig 2 sketches the overall prevalence in both patches with three different scenarios of residence 
times matrix P, both the constant and state-dependent case. The disease persists since the 
overall TZq > 1 in all three cases. 

3. The case where there is no movement between patches, that is, p \2 = P 21 = 0 ( pn = P 22 = 1) 
and 072 = <t 21 = 0 (or pi 2 (Ii, h ) = P 2 i(Ih h) = 0), corresponds to the case where the system 
behaves as two isolated patches. In this case the disease dies out or persists in Patch i if TZ l 0 is 
above or below unity in both approaches. This is illustrated on Fig 1(a) and Fig 1(b) where 
the disease dies out in Patch 1 ( Fig 1(a), blue solid line) where 7Zq = 3/777 = 0-^636 and 
the disease persists in Patch 2 (Fig 1(b), blue solid curve) where 7Zq = 37777 = 3.4286. For 
the state dependent case ( dashed blue curves on on Fig 1(a) and Fig 1(b)) the outcome is 
similar to the constant residence times case. 

4. In Fig 4(a) and 4(b), we explore the cases where there is symmetry (073 = Gji) with 
073 = Pij( 0,0). We supposed in this case that Patch 2 has higher risk (/?2 = 1.2) and Patch 
1 has lower risk {f3\ = 0.3). As can be intuitively deduced, the prevalence in Patch 1 is at 
its highest in the case of “high mobility” (072 = 021 = 1 ), and decreasing as 07 ,- decreases 
(with i 7 ! j). Conversely, prevalence in Patch 2 is at its highest under very “low mobility” 
(cr 12 = 021 = 0) and decreases as 077 increases. Note that 07 /, with j / j is proportional to 
Pij(I\ . I2) which is the actual residence time. 

5. We continue to explore the asymmetric case (073 / 037 ), that is, there is more mobility to¬ 
wards one patch. In Fig 5(a), the prevalence in Patch 1 (low risk) is at its highest if there is 
“high mobility” from Patch 1 to Patch 2 (072 = 1) and no mobility from Patch 2 to Patch 1 
(o '21 = 0 ), the prevalence decreases along with 072 . If the programmed response of residents 
of Patch 1 is to reduce their mobility (072 = 0) then, even if the mobility of residents in 
the high risk Patch 2 is extremely high (071 = 1), still the prevalence in Patch 1 is at its 
lowest. Similar remarks hold for Fig 5(b) regarding the prevalence in Patch 2 (high risk) 
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Figure 3. Dynamics of I\ and I2 for varying Oij for the state-dependent I2) where IZo < 1 . This is obtained 

by using the values f 3 \ = 0.2 and 02 = 0.3. In all the three cases, the disease dies out in both patches. The black 
curves represent the symmetric case ( p\2 = P21 = 0-12 = 021 = 0.5 ), the blue line represent the polar case ( 
P12 = P21 = o' 12 = 021 = 0) and red curves represent high mobility case (pi2 = p 21 = 0A2 = 021 = 1 )• 


under different mobility schemes. 

6 . Finally, Figure 6 presents the dynamics of the infected in both patches for the (conventional) 
case where p \2 = 0 ( and pu = 1). This case is particularly interesting since the residence 
time matrix P is not irreducible (hence the hypothesis of Theorem 2.1 fails) but 7 £q(P) = 
1.8929 > 1. As predicted by Theorem 2.2, the disease in Patch 2 is persistent. Also, it worth 
noticing that in Fig 6, I\ persists as well even though 7 £q(P) = 0.4455 < 1, as the condition 
7£g(P) > 1, for i = 1, 2, is sufficient but not necessary for persistence in Patch i. 


4. Final epidemic size 

Although the disease dynamics described here are not those of a controlled epidemiological system 
(the TZq is that corresponding to an uncontrolled system) they are still of considerable interest. 
The study of the role of residence time matrices on the dynamics of a single outbreak within a 
Susceptible-Infected-Recovered (with immunity) or SIR model without births and deaths is relevant 
to the development of public disease management measures [14, 26, 33]. Under the parameters and 
definitions introduced earlier, and making use of the same notation, we arrive at the following 
system of nonlinear differential equations: 


Si = - 


PiPu 


+ 


Pjpiij 


PnNi+PjiNj Pij N i+PjjNj 


Sih- 


fiiPiiPji 


h = 


PiPu 


PjPli 


PjPijPjj 


PiiNi+PjiNj Pij N i+Pjj N j 


hWji _I_I e. I. I I OiViiVji 1 PjPijPjj l Q T ■ _ a • T 

PiiNi+PjiNj PijNi+PjjNj J ^ \pnNi+ P jiNj ^ PijNi+PjjNj) 


Silj, 


(9) 
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(a) The level of prevalence in Patch 1 (low risk) (b) The level of prevalence in Patch 2 (high risk) 
seems to decrease as 012 and 021 decrease. seems to increase as (T 12 and 021 decrease. 


Figure 4 . Dynamics of 1 1 and I2 for varying aij for the state-dependent pij(Ii, I2) approach. 



to decrease as 012 and <721 decrease. seems to increase as 012 and 021 decrease. 


Figure 5 . Dynamics of and Ii and I2 for varying aij, but non-symmetric, for the state-dependent Pij(Ii, I2)■ 


where Ri denotes the population of recovered immune individuals in Patch i , cti is the recovery rate 
in Patch i and iVj = Si + R + Ri, for i = 1,2. 


The basic reproduction number TZq, is by definition the largest eigenvalue of 2 x 2 (n x n for the 
general case) next generation matrix, 


- FV _1 


01 P 11 


@2 p 1 2 ) 

1 Ni 

Pll A’i+p 2 lA r 2 


P12N1+P22N2 J 

] Ot 1 

P1P11P21 


P2P12P22 ] 

1 N2 

Pll M1+P21JV2 


Pl 2 -^1 ~\~P22 N2 ) 

1 Oil 


P1P11P21 , P2P12P22 \ Ni_ \ 

P11N1+P21N2 ' P12N1+P22N2 J a 2 I 

_ S1P21 j @ 2 P 22 j AU I 

Pll Vi+p 2 lA r 2 ' P12N1+P22N2 J a 2 ) 


It has been shown (see [34], for example) that not everybody gets infected during an outbreak, 
and so, estimating the size of the recovered population (the final epidemic size in the absence of 
deaths or departures) is tied in the solutions of the final size relationship, given in this case, by the 
system: 


log 

log 
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Kix 
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CN 
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( 10 ) 
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Figure 6 . Dynamics of h and I2 where pi 2 = 0 .In this case the residence time matrix P is not irreducible, the 
disease in Patch 2 persists nonetheless as predicted by the theorem 2 . 2 . 


where 


01 P 11 

+ 

faPl2 \ 

1 Ni 

P11-Y1+P21-/V2 

P12N1+P22 y 2 j 


ftlPllP21 

+ 

P2P12P22 ^ 

| A 5 
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P12-Y1+P22 N2 J 
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{ P1P21 _I_ 02P22__\ Y 2 

\Pll-Yi+p 2 lY 2 P12Y1+P22Y2 J 02 


The relationship (10) is obtained by using the fact that, in (9), we have Si + Ii = —Q-di > 0. 
This implies that lirn^oo /,(£) = 0 (for i = 1,2), since Si and I; are positive and integrating s l in 
(9), we obtain, after some tedious algebra Expression (10). The references [10, ] give more details 

on the computation of the final size relationship. 


It is important to observe that the next generation matrix and the matrix K defining the final 
epidemic size have the same eigenvalues. And so, the dominant eigenvalue, for both is TZo (although 
we note that we would not expect this to be the case in a controlled epidemiological system). 

The residence time matrix P plays an important role as evidenced by the dependence of the final 
epidemic size relation as in Fig 7. As we can notice in Fig 7, the prevalence in low risk Patch 1 is 
highest in the high mobility scheme where as in high risk Patch 2, the high mobility leads to the 
lowest prevalence. Also, as stated before ( lim Idt) = 0, for i = 1, 2.) with any typical outbreak 

t —>+00 

model, the disease ultimately dies out from both patches [35]. 


5. Conclusion and Discussions 

Heterogeneous mixing in multi-group epidemic models is most often defined in terms of group 
specific susceptibility and average contact rates captured multiplicatively by the transmission pa- 
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Prevalence in Patch 1 


Prevalence in Patch 2 




Figure 7 . The prevalence in patch 1 (low risk) reaches its highest when in extreme mobility case (solid blue line) 
and is lowest when there is no mobility between the patches. The opposite of this scenario happens in patch 2 (high 
risk). 


rameter ( 3 . However, contact rates, in general, cannot be measured in satisfactory ways for diseases 
like influenza, measles or tuberculosis, due to the difficulty of assessing the average number of con¬ 
tacts per unit of time of susceptible populations in different locations for varied activities. In this 
paper we propose the use of residence times in heterogeneous environments, as a proxy for “effective” 
contacts over a certain time window; and develop a multi-group epidemic framework via virtual 
dispersal where the risk of infection is a function of the residence time and local environmental risk. 
This novel approach eliminates the need to define and measure contact rates that are used in the 
traditional multi-group epidemic models with heterogeneous mixing. 

Under the proposed framework, we formulate a general multi-patch SIS epidemic model with 
residence times. We calculate the basic reproduction number TZq which is a function of a patch 
residence-times matrix P. Our global analysis shows that the model is robust in the sense that 
the disease dynamics depend exclusively on the basic reproductive number when the residence 
times matrix P is “constant” (Theorem 2.1). We proved that the disease free equilibrium is globally 
asymptotically stable (GAS) if the basic reproduction number 1Z$ < 1 and that a unique interior 
endemic GAS equilibrium exists if TZ o > 1. This results holds as long as the residence time matrix 
P is irreducible, that is, the graph of the patches is strongly connected. 

Our further analysis (Theorem 2.2) provide easily accessible insights on the impact of the res¬ 
idence matrix P on the levels of infection within each patch. Our results imply that the infection 
risk (measured by B) and the residence time matrix (P) can play an important role in the endemic 
at the patch level. More specifically, the right combinations of the environmental risk level (B) and 
dispersal behavior (P) can either promote or suppress infection for particular patches. This work 
complements the results of Theorem 2.1 regarding the robust dynamics under the assumption that 
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P is strongly connected, i.e., irreducible. For example, when Theorem 2.2 is applied to the two patch 
case, residents of Patch 1 visit Patch 2 but not conversely. 

These significant differences that emerges from the study of residence times models (P a “con¬ 
stant”) includes the possibility of studying disease dynamics in non strongly connected pacha con¬ 
figuration. In particular, we found conditions that allow us to characterize the patch-specific disease 
dynamics as a function of the time spend by residents and visitors to the patch of interest. This 
approach allowed us to classify patches as sources or sinks of infection, a role that depends on risk 
(£>) and mobility (P). 

We also explored the case where the entries of residence times matrix P are no longer constant 
but rather prevalence dependent. We noticed that whenever the residence times are negatively cor¬ 
related with risk then prevalence will be higher in the riskier patch but much lower than if the 
residence times were independent of health status. We ran carefully designed simulations to gain 
insights on the use of phenomenological modeling approach (System ( 8 )), since the mathematical 
analysis would be in general challenging. 

Our proposed framework has been applied to the context of a two-patch SIR single outbreak 
model with the concept of residence times. We derived the final epidemic size relationship in order 
to capture the size of the outbreak. Our results show that the residence time matrix P plays an 
important role which evidenced by the dependence of the final epidemic size relation as in Fig 7. 

In both conventional and phenomenological approaches to residence times used in this paper, 
humans behavior and responses to disease risk are automatic: P is constant and predefined functions 
of health status. Recent studies [30, 38, 39, 40, 51] have incorporated behavior as a feedback response 
coupled with the dynamics of the disease. A model of the decision to spend time in patch i = 1,2 
based on individuals’ utility functions that include the possibility of adapting to changing contagion 
dynamics in the above two patch setting, using previous work [30, 19], is the subject of a separate 
study. 
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Appendix A: Computation of TZo 

Proof. The general SIS model with residence time is described by the system ( 6 ) 

I = diag(!V — J)Pdiag(R)diag(lV)^ 1 P t / — diag(d/ + 7 j)I. 

The right hand member of the above system be can clearly decomposed as T + V where 


T = diag(A — I)Pdiag(£>)diag(IV) 1 P < 7 and V = —diag(d/ + 7 j)I 
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The jacobian at the DFE of F and V are giving by: 


F = DT 


= diag(iV)Pdiag(£>)diag(iV) 1 P t and V = V 


DFE 


= —diag(d/ + 7 /) 


DFE 


The basic reproduction number TZq is given by the spectral radius of the next generation matrix 
— FV^ 1 [27, 55]. Hence, we deduce that 

IZq = p{— diag(IV)Pdiag(0)diag(.ZV') - 1 P t V r_1 ) 

For the two-patch SIS model ( 6 ), IZq is the largest eigenvalue of the following matrix: 
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Appendix B: Proof of Theorem 2.1 


The proof uses the method in [ ] which is based on Hirsch’s theorem [37]. 
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Theorem B.l (Hirsch [37]). 

Let x = F(x) be a cooperative differential equation for which is invariant , the origin is an 
equilibrium, each DF(x ) is irreducible, and that all orbits are bounded. Suppose that 

x>y => DF(x) < DF(y) for all x,y. 

Then all orbits in K" tend to zero or there is a unique equilibrium p* in the interior o/K” and all 
orbits in M" tend to p*. 

Proof of Theorem 2.1. 

Equation ( 6 ) can be written as: 

I = (F + V)I — diag(/)Pdiag(S)diag(iV)“ 1 P t / (11) 

where F = diag(lV)Pdiag(B)diag(lV) _ 1 P t and V = —diag (dj + 7 /), as defined in Appendix A. Let 
us denote by X(I) the semi flow induced by (11). Hence 

DX(I) = diag(lV - 7)Pdiag(H)diag(lV)- 1 P t + V - W(h,I 2 ) (12) 

where W(I\ , I 2 ) = diag(Pdiag(H)diag(lV) _ 1 P*/). Since P is irreducible and I < N, DX(I) is clearly 
Metzler irreducible matrix. That means, the flow is strongly monotone. Plus, DX(I) is clearly 
decreasing with respect of I. Hence, by Hirsch’s theorem either all trajectories go to zero or go 
to an equilibrium point / ^> 0. From the relation (12), we have DX( 0) = F + V where F and V 
are the one defined previously in Appendix A. However, since F a nonnegative matrix and V is 
Metzler, we have the following equivalence 

a(F + V) < 0 <=► p(—FV~ l ) < 1 

where a(F + V) is the stability modulus, i.e: the largest real part of eigenvalues, of F + V 
and p(y—FV~ l ) the spectral radius of —FV -1 . Hence, the DFE is globally asymptotically stable 
if IZo = p(—FV~ X ) < 1. And if 1Z 0 > 1, i.e: a(F + V) > 0, the DFE is unstable [ ]. Since, we 
have proved that DX(I) is a Metzler matrix, to prove the local stability of the endemic equilibrium 
7 > 0, we only need to prove that it exists w 3> 0 such that DX(I)w < 0 [ ]. The endemic 
equilibrium I 3> 0 satisfies the equation 

(F + V)I - diag(7)Pdiag(H)diag(iV) _ 1 P t / = 0 

Hence, 

DX(I)I = -W(I)I < 0 

Hence, with w = I, we deduce that I is locally stable. With the attractivity of / guaranteed 
Hirsh’s theorem, we conclude that the endemic equilibrium I 0 is globally asymptotically stable 
if IZo > 1 . 

Finally, if IZo = 1, we have a(F + V) = 0. It exists c 2 > 0 such that (F + F) f c = 0. By 
considering the Lyapunov function V = ( c\I ). This function is definite positive and its derivation 
along the trajectories if ( 11 ) is 

L - (cl/) 

= (c\(F + V)I - diag(/)Pdiag(H)diag(lV) _ 1 P < /^ 

= — ^c|diag(/)Pdiag(H)diag(lV)"' 1 P f /^ 


< 


0 


(13) 
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Plus V = 0 only at the DFE. Hence the DFE is GAS if TZo = 1. This completes the proof of the 
theorem 2.1. □ 

Appendix C: Proof of Theorem 2.2 

Proof. Since Model (6) has the compact global attractor Cl, then according to Theorem (2.1), we 
can expect that lim^oo I,ft) < thus for time large enough, we can have ^ — I* > 0, therefore 
we have 






(di + '7 i)Ii 


which indicates follows when 72 q(P) > 1 


Then apply the average Lyapunov Theorem [45], we can conclude that liminf^oo I%(t) > 0, i.e., 
the disease in the residence Patch i is persistent if 7 Lq(P) > 1 . 


If pij > 0 and pkj = 0 for all k = 1,.., n, and k ^ i, this implies that if there is a portion of the 
residence Patch i population flowing into the residence Patch j, then there is no other residence 
Patch k where k j, i.e., 


which also implies that 
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then we can conclude that Model (6) can have an equilibrium since under these conditions, 
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Therefore, Ii = 0 is the invariant manifold for Model (6). 

On the other hand, when these conditions hold, then we have 




Pij d] 




= «jxE(f h>«- 
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Therefore, if 77 q(IP) = -Rq x i \jr) Pv < 1- then we have the following inequality: 


% = e"=i 




V n T7i • -fc- 
2^fc=l Pkl du 


(di T 


< Ii 
= Ii 


6 » 

di 2-‘j=1 


Pjpii 


Efe = 1 Pfei 


E?=i - (di + 7i) 


- (di + 7*) 


< 0 


Therefore, we have limt-*^ Jj(f) = 0, i.e., there is no endemic in the residence Patch i. 


□ 
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